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Abstract.  A  method  based  on  the  simple  shear  flow  modified  by  the  introduction  of  a  deterministic  non-conservative  force 
and  a  stochastic  process  is  proposed  to  measure  the  Navier-Stokes  shear  viscosity  in  a  granular  fluid  described  by  the  Enskog 
equation.  The  method  is  implemented  in  DSMC  simulations  for  a  wide  range  of  values  of  dissipation  and  density.  It  is  observed 
that,  after  a  certain  transient  period,  the  system  reaches  a  hydrodynamic  stage  which  tends  to  the  Navier-Stokes  regime  for 
long  times.  The  results  are  compared  with  theoretical  predictions  obtained  from  the  Chapman-Enskog  method  in  the  leading 
Sonine  approximation,  showing  quite  a  good  agreement,  even  for  strong  dissipation. 


1.  INTRODUCTION 

Under  the  assumption  of  molecular  chaos,  the  velocity  distribution  function  of  a  dense  fluid  of  inelastic  hard  spheres 
obeys  the  Enskog  kinetic  equation,  suitably  modified  to  account  for  the  collisional  inelasticity  through  a  constant 
coefficient  of  normal  restitution  a  [1].  By  an  extension  of  the  Chapman-Enskog  method,  the  Navier-Stokes  transport 
coefficients  for  a  single  fluid  have  been  derived  in  the  first  Sonine  approximation  as  nonlinear  functions  of  the  packing 
fraction  <j>  and  the  coefficient  of  restitution  [2].  The  shear  viscosity  for  a  dense  binary  mixture  has  also  been  evaluated 
within  the  same  Sonine  approximation  [3]. 

Given  that  the  first  Sonine  approximation  is  known  to  be  rather  accurate  in  the  case  of  elastic  hard  spheres  [4],  a 
natural  question  is  whether  this  approximation  is  still  reliable  for  finite  inelasticity.  This  point  is  of  interest  since  in 
the  Sonine  approximation  the  distribution  function  is  expanded  around  the  local  equilibrium  distribution.  However, 
a  granular  system  is  never  in  equilibrium  [5].  In  the  undriven  case,  the  uniform  reference  state  corresponds  to  the 
so-called  homogeneous  cooling  state  (HCS),  in  which  the  time  dependence  only  occurs  through  the  temperature.  The 
velocity  distribution  function  of  the  HCS  deviates  from  a  Gaussian,  as  measured  by  a  non-zero  value  of  the  fourth- 
cumulant  and  by  an  overpopulated  high-energy  tail  [6].  This  fact  could  cast  some  doubts  about  the  accuracy  of  the 
Sonine  approximation  to  compute  the  transport  coefficients. 

The  shear  viscosity  T]  is  perhaps  the  most  widely  studied  transport  coefficient  in  a  granular  fluid.  For  a  low-density 
gas  ((/)  — >  0),  Brey  et  al.  [7]  measured  r\  from  DSMC  by  analyzing  the  time  decay  of  a  weak  transverse  shear  wave, 
and  observed  a  qualitative  good  agreement  with  the  Sonine  prediction.  However,  while  this  method  is  physically  quite 
elegant,  it  is  not  perhaps  sufficiently  efficient  to  measure  rj  accurately.  More  recently,  the  low-density  Navier-Stokes 
transport  coefficients  of  the  HCS  have  been  measured  from  DSMC  by  using  the  associated  Green-Kubo  formulas 
[8].  In  the  case  of  a  system  heated  by  the  action  of  an  external  driving  or  thermostat,  the  associated  shear  viscosity 
(which  is  different  from  that  of  the  HCS)  has  been  computed  from  the  Chapman-Enskog  method  in  the  first  Sonine 
approximation  and  from  DSMC,  for  a  dilute  single  gas  [9]  as  well  as  for  a  dense  binary  mixture  [3,  10].  On  the  other 
hand,  to  the  best  of  our  knowledge,  the  shear  viscosity  of  the  HCS  for  a  dense  gas  described  by  the  Enskog  equation 
has  not  been  measured  from  DSMC  before.  This  issue  is  the  main  goal  of  this  paper.  As  a  second  point,  we  want  to 
address  the  question  of  whether  a  granular  fluid  reaches  a  hydrodynamic  regime  that  can  be  described  by  the  Navier- 
Stokes  equations,  even  far  from  the  quasi-elastic  limit.  In  fact,  a  granular  fluid  in  an  inhomogeneous  steady  state  is 
inherently  non-Newtonian  [5,  1 1],  so  that  a  Newtonian  regime  requires  unsteady  states  sufficiently  separated  from  the 
initial  preparation. 

In  order  to  investigate  the  problem,  we  have  carried  out  numerical  simulations  of  the  Enskog  equation  by  means  of 
an  extension  of  the  DSMC  method  [12]  for  0.6  <  a  <  1  and  0  <(j)  <  0.5  and  have  considered  the  so-called  simple  shear 
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flow  state.  To  allow  the  granular  fluid  to  approach  a  Newtonian  regime,  an  external  “friction"  force  with  a  negative 
friction  coefficient  is  applied  to  compensate  for  the  inelastic  cooling,  so  that  viscous  heating  produces  a  monotonic 
decrease  of  the  Knudsen  number  (reduced  shear  rate)  with  time.  In  addition,  to  mimic  the  conditions  appearing  in 
the  Chapman— Enskog  method  to  Navier-Stokes  order,  a  stochastic  process  is  introduced,  according  to  which  every 
particle  has  a  certain  probability  per  unit  time  to  have  its  velocity  replaced  by  a  new  velocity  sampled  from  the  HCS 
distribution. 

Our  simulation  results  confirm  that,  regardless  of  its  initial  preparation,  the  fluid  reaches  after  a  few  collisions  per 
particle  a  hydrodynamic  state  whose  temporal  evolution  is  governed  by  that  of  the  granular  temperature.  Moreover, 
when  the  shear  viscosity  is  measured  as  the  ratio  between  the  shear  stress  and  the  shear  rate  for  long  times,  it  agrees 
reasonably  well  with  the  theoretical  expressions  derived  in  the  first  Sonine  approximation  [2],  the  deviations  being 
comparable  to  or  smaller  than  the  ones  in  the  elastic  case. 


2.  MODIFIED  SIMPLE  SHEAR  FLOW 

The  simple  shear  flow  is  macroscopically  characterized  by  a  constant  density  n,  a  uniform  granular  temperature  7  ,  and 
a  linear  velocity  field  u(r)  =  a  •  r,  where  the  rate  of  strain  tensor  a  is  «,/  =  a8jx§jy,  a  being  the  constant  shear  rate  [13]. 
At  the  level  of  the  one-body  distribution  function,  the  simple  shear  flow  becomes  uniform  in  the  Lagrangian  frame  of 
reference,  i.e.,  f(r. \.t)  =  f(\(r).t),  where  V(r)  =  v  —  u(r).  Consequently,  the  Enskog  equation  in  this  problem  reads 

dtf-aVy^-f  +  F{f}=J[f,f}.  (1) 

Here  J[f,f]  is  the  Enskog  operator,  which  in  the  simple  shear  flow  is  given  by 

J[f,f]  =  v2x(n)  jd\ !  j  da&(a  •  g)(cr  •g)[a_2/(V,)/(V,1  +  a  •  <y)  -  /(V)/(Vi  -  a  •  a)],  (2) 

where  a  is  the  diameter  of  a  sphere,  x{n)  is  the  contact  value  of  the  pair  distribution  function,  a  is  a  unit  vector  directed 
along  the  centers  of  the  two  colliding  spheres  at  contact,  g  =  V  —  Vi  is  the  relative  velocity,  a  <  1  is  the  coefficient 
of  normal  restitution,  a  =  ac .  and  the  primes  denote  pre-collisional  velocities,  namely  V'  =  V—  ^(1+  a  ~ 1 )  ( <T  •  g)£T 
and  Vj  =  Vi  +  |(1  +  «  1 ) (O'  •  g)0\  Finally,  F[f\  in  Eq.  (1)  is  an  operator  representing  some  type  of  external  action 
on  /,  absent  in  the  true  simple  shear  flow  problem,  that  will  be  chosen  later  on.  This  extra  term  is  assumed  to  verify 
the  conditions 

J  dV{l,\,V2}F\f\  =  {0,0,-3nYT/m},  (3) 

so  that  it  does  not  affect  the  mass  and  momentum  conservation  laws,  but  in  general  contributes  to  the  energy  balance 
equation  through  the  “heating”  rate  y.  The  time  evolution  for  the  granular  temperature  is 

dtT  =  -^aPxy-CT  +  yT,  (4) 

where  Pxy  =  Pxy  +  P(v  is  the  total  shear  stress  and  C  is  the  cooling  rate, 

pfj  =  J  dWmV,Vjf(y ),  Pfj  =  ]-^ma3x  j  doofr  J  d\  J  </Vi0(ff  •  g)(ff -g)2/(V  +  a  •  (5) 

2 

C=^(1-«2)z  JdoJdyJ  dV  i®(o  •  g)(cr  •  g)3/(V  +  a  •  <y)/(Vi).  (6) 

The  first  term  on  the  right-hand  side  of  Eq.  (4)  represents  viscous  heating  effects,  the  second  term  corresponds  to  the 
cooling  due  to  the  inelasticity  of  collisions,  while  the  third  term  is  the  contribution  due  to  the  external  action  F\f\. 

In  the  true  simple  shear  flow,  i.e.,  with  F[f  \  =  0,  the  temperature  evolution  is  governed  by  the  competition  between 
the  viscous  heating  and  the  collisional  cooling  effects  until  a  steady  state  is  reached  when  both  effects  cancel  each 
other.  However,  this  steady  state  is  inherently  non-Newtonian  [1 1],  so  that  the  Navier-Stokes  shear  viscosity  coefficient 
cannot  be  measured,  even  in  the  quasi-elastic  limit.  The  domain  of  validity  of  the  Newtonian  description  is  restricted 
to  small  values  of  the  Knudsen  number  Kn  =  X / 1,  where  A  =  \ / \/2TCncj2 x(n)  is  the  mean  free  path  and  f  is  the 
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characteristic  hydrodynamic  length.  In  the  present  problem,  the  only  hydrodynamic  length  is  t  =  \JTT /m /a.  In  the 
steady  shear  flow  Kn  is  a  function  of  a  that  cannot  be  controlled  independently.  In  order  to  reach  asymptotically  small 
values  of  Kn  for  any  value  of  a  we  need  to  avoid  a  steady  state  and  have  a  monotonically  increasing  temperature. 
Consequently,  we  must  modify  the  original  shear  flow  problem  by  introducing  an  external  driving  mechanism, 
represented  by  F\f\,  which  exactly  compensates  for  the  collisional  cooling  term  in  Eq.  (4).  Specifically,  the  heating 
rate  introduced  in  Eq.  (3)  is  chosen  as  y  =  £,  so  dtT  =  —(2/3 n)aPxy. 

So  far,  apart  from  the  requirement  y  =  £,  the  explicit  form  of  F\f\  remains  still  open.  Here  we  will  determine 
F\f\  by  requiring  that  the  kinetic  equation  (1)  to  first  order  in  Kn  be  the  same  as  the  one  found  by  applying  the 
Chapman-Enskog  method  for  states  close  to  the  (local)  HCS  [2],  To  that  end  we  formally  expand  f  as 

/(V)=/(0)(V)+/(1)(V)  +  ---,  (7) 

where  ftkl  (V)  is  of  order  k  in  Kn.  Moreover,  the  time  dependence  of  f  occurs  entirely  through  the  temperature.  Note 
that,  by  definition,  a  ~  Kn.  On  physical  grounds,  Pxy  is  at  least  of  first  order  in  Kn.  Therefore,  given  that  y  =  £, 
dtT  ~  Kn2  and  so  d,f  can  be  neglected  to  first  order.  Inserting  Eq.  (7)  into  Eq.  (1),  we  get,  through  first  order  in  Kn. 

F[/(0)]=/0)[/(0),/(0)],  (8) 

F[/«]  -aVyJ^fW  =  -L[f^]+aAy  [df^/dVx]  ,  (9) 

where 

j^[x1Y]  =  a2x{n)  f  d\\  J do®(d-%)(o-%)[a-2x(y')Y(y\)-x(y)Y(yy)},  do) 

L[X]  =  -7^[/0\X]-/°)[X,/0)],  (11) 

Ai[X]=e3x(n)  j dVi  j  Ja©(a-g)(a-g)a,[a-2/(0)(V,)^(V,1)+/(0)(V)Z(V1)].  (12) 

We  require  that  fl0>  be  the  HCS.  This  implies  that  [6] 

^[/(0)]  =  [v/(0)(V)]  ,  (13) 

where  Cl{})  is  the  cooling  rate  in  the  HCS,  which  is  obtained  from  Eq.  (6)  by  setting  a  =  0  and  /  — >  f-y> .  Next,  we 
consider  the  linear  integral  equation  (9)  for  /* 1  f .  This  equation  coincides  with  the  one  derived  in  Ref.  [2]  from  the 
general  Chapman-Enskog  method  specialized  to  Vn  =  VT  =  V  •  u  =  0,  provided  that 

F[fW]  =  ~^0)TdTfW  =  [V/(1)(V)]  +  ^C(0)/(1)(V),  (14) 

where  in  the  last  step  we  have  taken  into  account  that  Kn  °c  T  12  and,  by  dimensional  analysis,  /'" 1  •'  ( V)  = 
«Vq  3<I>(V/vo)Kn,  where  vo  =  y/2T /m  is  the  thermal  velocity,  <t>  being  a  dimensionless  function  of  the  reduced  ve¬ 
locity.  The  simplest  choice  for  F\f\  compatible  with  Eq.  (3)  (with  y  =  £)  and  with  Eqs.  (13)  and  (14)  is 

F[/]  =  ^^-(V/)  +  ic(/-/(0)).  (15) 

In  summary,  our  modified  simple  shear  flow  problem  is  described  by  Eq.  (1)  along  with  the  choice  (15).  The  first 
term  on  the  right-hand  side  of  Eq.  (15)  represents  the  effect  of  a  deterministic  non-conservative  force  of  the  form 
2  £mV,  which  does  work  to  compensate  for  the  collisional  energy  loss.  The  shear  viscosity  of  a  granular  fluid  mixture 
when  only  this  term  is  present  in  F\f\  has  been  determined  from  the  Chapman-Enskog  method  and  measured  in 
DSMC  simulations  [3,  10].  However,  this  term  does  not  suffice  to  get  the  Navier— Stokes  shear  viscosity  of  the  HCS, 
but  it  must  be  supplemented  by  the  second  term  on  the  right-hand  side  of  Eq.  (15).  The  latter  term  represents  the  action 
of  a  stochastic  process,  whereby  each  particle  has  a  probability  per  unit  time  equal  to  /  £  of  changing  its  velocity  by 
a  new  velocity  sampled  from  the  (local)  HCS  distribution  /(°) .  When  this  stochastic  term  is  moved  to  the  right-hand 
side  of  the  Enskog  equation  (1),  it  can  be  interpreted  as  a  BGK-like  relaxation  term. 
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FIGURE  1.  Left  panel:  Parametric  plot  of  the  reduced  normal  stresses  Pj’j/nT  versus  Kn2  for  a  =  0.8  and  <p  =  0.5.  Right  panel: 
Parametric  plot  of  the  fourth  cumulant  c,  relative  to  that  of  the  HCS,  versus  Kn2  for  a  =  0.9  and  <j)  =  0.2.  0.3,  and  0.4. 


3.  DSMC  RESULTS 


We  have  solved  the  Enskog  equation  (1)  with  the  presence  of  F\f\  as  given  by  Eq.  (15)  by  means  of  an  extension  to  the 
Enskog  equation  [12]  of  Bird’s  DSMC  method.  Since  the  system  is  uniform  in  the  Lagrangian  frame,  there  is  no  need 
to  split  it  into  cells.  In  the  implementation  of  F[f]  the  cooling  rate  is  self-consistently  measured  in  the  simulations. 
The  stochastic  term  is  simulated  by  randomly  choosing  during  a  small  time  step  St  a  fraction  \C,8t  of  particles;  the 
original  velocity  V0id  of  each  one  of  those  particles  is  replaced  by  a  new  velocity  Vnew  =  wT  C,  where  C  is  the  velocity 
of  a  particle  in  a  reservoir  kept  at  the  HCS  normalized  to  unit  temperature. 

The  shear  rate  and  the  initial  temperature  are  such  that  Kn  =  0.1  at  t  =  0  for  all  density  and  coefficient  of  restitution. 
Besides,  the  initial  velocity  distribution  is  a  Maxwellian.  In  the  course  of  the  simulations,  the  kinetic  and  collisional 
transfer  contributions  (5)  to  the  pressure  tensor  are  evaluated.  From  the  shear  stress  Pxy(t)  the  shear  viscosity  is 
measured  as  a  function  of  time  as  T){t)  =  —Pxv(t)/a.  Since  the  Knudsen  number  Kn(f)  °c  1  /  y/Tjt)  monotonically 
decreases  with  increasing  time,  the  Navier-Stokes  shear  viscosity  is  identified  as  rj  (r)  for  long  times. 

As  said  in  Section  1,  our  main  objective  is  to  compare  the  kinetic  (r)k)  and  total  (77)  shear  viscosity  measured  in  the 
simulations  with  the  expressions  derived  from  the  Chapman-Enskog  method  by  using  the  first  Sonine  approximation 
f{  1 1  — >  — a(mr\k / nT2) fMVxVy ,  /+  being  the  Maxwellian  distribution.  The  expressions  are  [2] 


rik(oc,  0)  =  t]o 


_ l-g^X(^)(l  +  «)(l-3«) _ 

354XOXI  + a)  [16(13  — a)  -  3(4  -  3a)c0(a)]  ’ 


(16) 


n(«0)  =  n  («0) 


l  +  ^OxOXl  +  a) 


384  9  , 

+  T?0  ^O^OXi  +  a) 


1  -  32C°(a) 


(17) 


where  T]o  =  ( 5 / 1 6 rr 2 )  y7/ nT ]  n  is  the  shear  viscosity  of  a  dilute  gas  in  the  elastic  limit,  (j)  =  |n(73  is  the  packing 
fraction,  and  cq  is  the  fourth  cumulant  of  the  HCS.  A  good  estimate  of  cq  is  [6] 


c0(a) 


32(1  — a)(l  — 2a2) 
81-  17a  +  30a2(l -a)' 


(18) 


In  addition,  we  take  the  Carnahan— Starling  approximation  x(<j>)  =  (1  —  0 / 2) /( 1  —  <j>)3. 

Figure  1  shows  the  kinetic  part  of  the  normal  stresses  Pjj  (f),  relative  to  nT  (t),  and  the  fourth  cumulant  c(t),  relative  to 
its  HCS  value  cq,  as  functions  of  the  time-dependent  Knudsen  number  Kn(f)  =  A / f  (T)  =  act / \2<j> xi<f>)  T {t) / in,  for 
some  representative  cases.  Note  that,  as  time  grows,  the  Knudsen  number  Kn  monotonically  decreases  from  its  initial 
value  Kn  =  0.1,  behaving  as  Knit)  ~  t  1  for  asymptotically  long  times.  We  observe  that  in  the  long-time  limit,  i.e., 
for  Kn  — >  0,  the  system  tends  to  an  isotropic  state  with  a  fourth  cumulant  c  equal  to  that  of  the  HCS.  This  supports  the 
expectation  that  the  asymptotic  state  of  our  modified  simple  shear  flow  is  the  HCS,  despite  the  fact  that  the  temperature 
is  increasing  rather  than  decreasing,  in  agreement  with  Eq.  (8).  It  is  worth  remarking  that,  according  to  Fig.  1,  after  a 
short  transient  period  the  fluid  reaches  a  hydrodynamic  regime  where  the  normal  stresses  and  the  cumulant  are  linear 
functions  of  Kn2  (Burnett-order  effects). 
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FIGURE  2.  Parametric  plot  of  the  reduced  kinetic  shear  viscosity  rj^(a,0)/t7^(l,0)  versus  Kn2  for  a  =  0.8  and  (j)  =  0.2.  The 
dotted  line  represents  the  estimated  Navier-Stokes  value. 


As  an  illustration  of  how  the  Navier-Stokes  shear  viscosity  is  evaluated  from  DSMC,  Fig.  2  shows  the  time- 
dependent  kinetic  shear  viscosity  T]k{t)  =  —Pky(t) /a,  relative  to  its  (time-dependent)  theoretical  Sonine  value  in 
the  elastic  limit,  as  a  function  of  Kn2(f)  for  the  case  a  =  0.8,  <p  =  0.2.  Figure  2  clearly  shows  that  the  ratio 
T]k(a,(j))/rik(l,(j))  reaches  a  plateau  for  long  times  (small  Knudsen  numbers)  that  can  be  identified  as  the  Navier- 
Stokes  value.  The  same  procedure  has  been  followed  to  measure  the  kinetic  and  collisional  parts  of  the  Navier-Stokes 
shear  viscosity  for  different  values  of  dissipation  and  density. 


FIGURE  3.  Plots  of  the  kinetic  part  of  the  shear  viscosity  (left  panel)  and  of  the  total  shear  viscosity  (right  panel),  relative  to  their 
respective  theoretical  values  in  the  elastic  limit  in  the  first  Sonine  approximation,  as  functions  of  the  packing  fraction  for  a  =  0.9, 
0.8,  and  0.6.  The  symbols  are  simulation  results,  while  the  lines  are  the  theoretical  predictions  in  the  first  Sonine  approximation. 


FIGURE  4.  Plots  of  the  kinetic  part  of  the  shear  viscosity  (left  panel)  and  of  the  total  shear  viscosity  (right  panel),  relative  to 
their  respective  theoretical  values  in  the  elastic  limit  in  the  first  Sonine  approximation,  as  functions  of  the  coefficient  of  restitution 
for  <p  =  0,  0.2,  and  0.4.  The  symbols  are  simulation  results,  while  the  lines  are  the  theoretical  predictions  in  the  first  Sonine 
approximation. 
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The  density  dependence  of  T]k  and  rj  for  three  values  of  the  coefficient  of  restitution  is  displayed  in  Fig.  3,  while 
Fig.  4  shows  the  influence  of  dissipation  for  three  values  of  the  packing  fraction.  We  observe  a  general  good  agreement 
between  DSMC  results  and  the  predictions  from  the  first  Sonine  approximation,  even  for  strong  dissipation  and  high 
densities.  Both  theory  and  simulation  show  that,  at  a  given  value  of  a,  r\k(a,ty)  >  r; *(1,0)  if  the  packing  fraction 
is  smaller  than  a  certain  value  0,* (a),  while  rjk(a,tj))  <  rj*(  1 , 0)  if  0  >  0q (a).  A  similar  behavior  occurs  for  the 
total  shear  viscosity  with  a  different  value  0o(a).  The  influence  of  a  on  both  </>/]  and  0, q  is  rather  weak;  according 
to  Eqs.  (16)— (18),  (0Q, 0o)  =  (0.12,0.09),  (0.13,0.09),  and  (0.15,0.10)  for  a  =  0.9,  0.8,  and  0.6,  respectively.  As 
a  consequence,  while  in  a  dilute  granular  gas  (0  <  0.1)  the  shear  viscosity  increases  with  inelasticity,  the  opposite 
happens  for  sufficiently  dense  fluids  (0  >  0. 1). 


4.  CONCLUDING  REMARKS 

In  this  paper  we  have  proposed  a  method  to  measure  the  Navier-Stokes  shear  viscosity  of  a  moderately  dense  granular 
fluid  described  by  the  Enskog  equation.  The  idea  is  to  consider  the  simple  shear  flow  (which  is  uniform  in  the 
Lagrangian  frame),  modified  by  the  presence  of  a  deterministic  non-conservative  force  (which  compensates  for  the 
collisional  cooling)  along  with  a  stochastic  BGK-like  term.  Under  these  conditions  the  Knudsen  number  of  the  problem 
decreases  with  increasing  time,  so  that  the  system  reaches  a  hydrodynamic  Navier-Stokes  regime  for  long  times.  This 
procedure  allows  one  to  evaluate  the  Navier-Stokes  shear  viscosity  in  an  efficient  way  by  means  of  the  DSMC  method. 
The  simulation  results  have  been  compared  with  predictions  from  the  Chapman-Enskog  expansion  in  the  first  Sonine 
approximation  [2].  The  results  show  that  the  Sonine  predictions  compare  quite  well  with  the  simulation  data  for  the 
wide  range  of  dissipation  (a  >  0.6)  and  density  (0  <  0.5)  explored.  This  agreement  is  significant,  given  that,  in  contrast 
to  the  elastic  case,  the  reference  state  is  not  a  (local)  Maxwellian  but  the  (local)  HCS,  which  exhibits  non-Maxwellian 
features,  such  as  a  non-zero  fourth  cumulant  and  an  overpopulated  high-velocity  tail  [6].  It  is  interesting  to  remark  that 
the  accuracy  of  the  Sonine  approximation  found  here  is  comparable  to  the  one  observed  for  elastic  fluids,  so  that  one 
could  expect  to  improve  the  agreement  by  considering  more  terms  in  the  Sonine  polynomial  expansion.  Finally,  we 
want  to  stress  that  the  system,  after  a  short  transient  period,  achieves  a  hydrodynamic  stage  prior  to  the  Navier-Stokes 
regime.  This  stage  could  be  used  to  study  nonlinear  transport  properties  (e.g.,  shear  thinning  and  viscometric  effects), 
although  this  issue  is  beyond  the  scope  of  this  paper. 
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